# load packages
library(ggplot2)
library(cregg)
library(foreign)
library(lmtest)
library(sandwich)
library(msm)
library(dplyr)
library(tidyverse)
library(hrbrthemes)
library(estimatr)
library(ggrepel)
library(stargazer)
library(here)
library(tidyr)
library(FindIt)

#set path to conjoint replication folder
setwd("C:/Users/ddow1/Dropbox/Community Policing MPP Paper/Code/replication/")

#Load data
e1 <- read_rds("mpp_cjt_data.Rds")

## Labels for plotting
cjt.labels = c(EconLoss = "ECONOMIC LOSSES", 
               Violence = "GENERAL VIOLENCE", 
               Effective = "EFFECT OF HELPING POLICE", 
               Neighbors = "NEIGHBORS HELP", 
               Reprisal = "CHANCES OF REPRISALS", 
               Police = "POLICE ORIGINS")

#Plot marginal means
ef1 <- chosen ~ Police + Reprisal + Violence + EconLoss + Neighbors + Effective
mms <- plot(mm(e1, ef1, id = ~instanceid, feature_labels = cjt.labels, 
            level_order = c("descending")), vline = 0.5, size = 3)


my_faces <- c(rep("plain", 2), "bold",
              rep("plain", 2), "bold",
              rep("plain", 2), "bold",
              rep("plain", 2), "bold",
              rep("plain", 3), "bold",
              rep("plain", 2), "bold")


pmm <- plot(mms, vline = 0.5, header_fmt = "%s", size = 2) + ggplot2::theme(
  legend.position = "none",
  panel.grid.major = ggplot2::element_blank(),
  panel.grid.minor = ggplot2::element_blank(),
  axis.text.y = element_text(face = my_faces, size = 11)) +
  ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower",
                                              xmax = "upper"),
                          size = 2, height = 0, na.rm = TRUE,
                          position = ggstance::position_dodgev(height = 1)) +
  ggtitle("Figure 1: Attribute Choice for Police Cooperation (Marginal Means)")

pmm

ggsave(filename = "fig_1_mm.pdf", 
       device = cairo_pdf, width = 8, height = 10)

#Plot AMCEs
amces <- cj(e1, ef1, id = ~instanceid, feature_labels = cjt.labels, 
  level_order = c("descending"))

amce.main <- plot(amces, size = 3)

pamce <- plot(amce.main, vline = 0.5, header_fmt = "%s", size = 2) + 
  ggplot2::theme(
  legend.position = "none",
  panel.grid.major = ggplot2::element_blank(),
  panel.grid.minor = ggplot2::element_blank(),
  axis.text.y = element_text(face = my_faces, size = 11)) +
  ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower",
                                              xmax = "upper"),
                          size = 2, height = 0, na.rm = TRUE,
                          position = ggstance::position_dodgev(height = 1)) +
  ggtitle("Figure 2: Average Marginal Component Effects (AMCE)")

pamce

ggsave(filename = "fig_2_amce.pdf", 
       device = cairo_pdf, width = 8, height = 10)


# HTE ANALYSIS for POLICE LEGITIMACY
#subgroup analysis for perceptions of police legitimacy

## 1) High = strong agree and agree vs 0) low = disagree, strong disagree
## drop neutral category, don't know, refuse
e1$pol_legit_lh <- e1$pol_legit
e1$pol_legit_lh <- na_if(e1$pol_legit_lh, '888')
e1$pol_legit_lh <- na_if(e1$pol_legit_lh, '999')
e1$pol_legit_lh <- na_if(e1$pol_legit_lh, '3')
e1$pol_legit_lh <- ifelse(e1$pol_legit_lh < 3, 1, ifelse(e1$pol_legit_lh > 3, 0, NA))
e1$pol_legit_lh <- factor(e1$pol_legit_lh, 0:1, c("Low Police Legitimacy", "High Legitimacy"))

## Marginal Means by Legitimacy Level
e1.legitsplitlh <- cj(e1, ef1, 
                      na.rm = T, id = ~ instanceid, estimate = "mm", h0 = 0.5, by = ~ pol_legit_lh, 
                      feature_labels = cjt.labels)
                        
plot(e1.legitsplitlh, feature_labels = cjt.labels, group = "pol_legit_lh", vline = 0.5, size = 2)

ggsave(("fig_a6_legit_mm.pdf"), 
       device = cairo_pdf, width = 8, height = 10)

## Contrasts in AMCEs by Legitimacy Level
legit = 
  e1 %>%
  dplyr::select(chosen, Police, Effective, Reprisal, 
         Violence, EconLoss, Neighbors, 
         instanceid, pol_legit_lh)

amce.legitlh <- cj(na.omit(legit), ef1, id = ~instanceid, estimate = "amce", 
                   by = ~pol_legit_lh, feature_labels = cjt.labels)
                  
diff_amce.legitlh <- cj(na.omit(legit), ef1, id = ~instanceid, estimate = "amce_diff", 
                        by = ~pol_legit_lh, feature_labels = cjt.labels)


plot(rbind(amce.legitlh, diff_amce.legitlh), size = 2) + ggplot2::facet_wrap(~BY, ncol = 3L) +  
  ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower", xmax = "upper"), 
                          size = 1, height = 0, na.rm = T) +
  ggtitle("Figure A7: AMCE by Perceptions of Police Legitimacy")

ggsave(filename = "fig_a7_amce_contrast.pdf", 
       device = cairo_pdf, width = 10, height = 8)


###############################################
##      Diagnostics on Conjoint
###############################################

task_anova = cj_anova(e1, ef1, by = ~cjt_task)

cleanup = tribble(~feature, ~label, 
                  "EconLoss", "Econ Loss", 
                  "Violence", "Level Violence", 
                  "Effective", "Pol. Effectiveness", 
                  "Neighbors", "Neighbor Cooperate", 
                  "Reprisal", "Risk of Reprisals",
                  "Police", "Police Origins")

stargazer(task_anova, type = "latex", summary = F, 
          title = "Table A3: Diagnostic: are respondent preferences changing across conjoint tasks? Nested model comparison test.", style = "ajps", 
          out = "table_a3_diagnostic.tex")

## Figure A2: Marginal Means by Conjoint Task
e1$cjt_task <- factor(e1$cjt_task)

task_prefs = cj(as.data.frame(e1), ef1,
                id = ~instanceid, by = ~cjt_task, estimate = "mm")

task_prefs %>% 
  left_join(cleanup) %>% 
  ggplot(aes(x = level, y = estimate, ymin = lower, ymax = upper, label = level, 
             shape = BY, color = BY)) + 
  geom_pointrange(position = position_dodge(width = .5), alpha = .8) + 
  coord_flip() +
  theme_light(base_family = "IBM Plex Sans Condensed") + 
  theme(panel.grid.major = element_blank(), legend.position = "top",
        strip.background =element_rect(fill = "black"),
        strip.text =element_text(color = "white", face = "bold", size = 4)) + 
  labs(y = "Estimate", 
       x = NULL, 
       color = "Task number:", 
       shape = "Task number:") + 
  ggtitle("Figure A2: Marginal Means by Conjoint Task") +
  scale_y_percent() + 
  scale_color_brewer(palette = "Dark2") +
  facet_grid(rows = vars(label), scales = "free", switch = "x")

ggsave(filename = "fig_a2.pdf",
       device = cairo_pdf, width = 5, height = 6)


#############################
### Causal Interactions
#############################

# Re-word neighbor values for clarity in plot
library(plyr)
e1$Neighbors <- revalue(e1$Neighbors, c("Will probably help"="Others probably help", "Will probably not help"="Others probably won't help"))
e1$Reprisal <- revalue(e1$Reprisal, c("No chance" = "No reprisal chance", 
                                      "Low chance"="Low reprisal chance", 
                                      "High chance"="High reprisal chance"))

e1$EconLoss <- factor(e1$EconLoss,ordered=FALSE,levels=c("Small","Large"))
e1$Violence <- factor(e1$Violence,ordered=FALSE, c("Low", "High"))
e1$Police <- factor(e1$Police,ordered=FALSE, c("From another place", "From here"))
e1$Neighbors <- factor(e1$Neighbors,ordered=FALSE, c("Others probably won't help","Others probably help"))
e1$Effective <- factor(e1$Effective,ordered=FALSE, c("Not help", "Help some"))
e1$Reprisal <- factor(e1$Reprisal,ordered=TRUE, c("No reprisal chance", "Low reprisal chance", "High reprisal chance"))


## Figure A3 Reprisal X Level of Violence
fit.i <- CausalANOVA(formula=ef1, int2.formula = ~ Reprisal:Violence,
                    data=e1, pair.id=e1$instanceid, diff=TRUE,
                    cluster=e1$instanceid, nway=2)


pdf('fig_a3_reprisal_conditional.pdf')
par(mgp=c(3,1.5,0))
plot(fit.i, type="ConditionalEffect", fac.name=c("Reprisal","Violence"), xlim = c(-0.1,0.1))
dev.off()


## Average Marginal Interaction Effects (AMIE)
e1$contestpair = 
  paste(e1$instanceid, e1$cjt_task)

# split sample into train and test data
set.seed(4600)

train.id = sample(unique(e1$instanceid), 
                  length(unique(e1$instanceid))/2, replace = F)
test.id = base::setdiff(x = unique(e1$instanceid), y = train.id)

train = e1[is.element(e1$instanceid, train.id), ]
test = e1[is.element(e1$instanceid, test.id), ]

fit <- CausalANOVA(chosen ~ EconLoss + Violence + Effective + Neighbors + Reprisal + Police, 
                   data=train, 
                   pair.id=train$contestpair,
                   diff=TRUE,
                   screen=TRUE, 
                   collapse=TRUE,
                   cluster=e1$instanceid, 
                   nway=2)
summary(fit)
## refit with test.CausalANOVA
fit.new <- test.CausalANOVA(fit, 
                            newdata=test, 
                            diff=TRUE,
                            pair.id=test$contestpair, 
                            cluster=test$instanceid)
summary(fit.new)

#Plot Figure A5
pdf('fig_a5_conditional.pdf', 
    title ="Figure A5")
par(mgp=c(3,1.5,0))
plot(fit.new, type="ConditionalEffect", 
     fac.name=c("Neighbors", "Police"))
dev.off()

ConditionalEffect(fit.new, treat.fac="Neighbors", cond.fac="Police")

#plot AMIEs
amies = rbind(fit.new$CI.table[[7]], 
              fit.new$CI.table[[8]], 
              fit.new$CI.table[[9]]) %>% 
  dplyr::mutate(label = rownames(.)) %>% 
  dplyr::filter(AMIE != 0)

amies$label

# clean up names
amies$label = stringr::str_replace_all(amies$label, 'EffectiveHelp some', 'Effect - help')
amies$label = stringr::str_replace_all(amies$label, 'EffectiveNot help', 'Effect - no help')
amies$label = stringr::str_replace_all(amies$label, "NeighborsOthers probably won't help", 'Neighbors - not help')
amies$label = stringr::str_replace_all(amies$label, 'NeighborsOthers probably help', 'Neighbors - help')
amies$label = stringr::str_replace_all(amies$label, 'ViolenceHigh', 'Violence - high')
amies$label = stringr::str_replace_all(amies$label, 'ViolenceLow', 'Violence - low')
amies$label = stringr::str_replace_all(amies$label, 'PoliceFrom another place', 'Police - from elsewhere')
amies$label = stringr::str_replace_all(amies$label, 'PoliceFrom here', 'Police - from here')
amies$label = stringr::str_replace_all(amies$label, 'EconLossLarge', 'Econ Losses - Large')
amies$label = stringr::str_replace_all(amies$label, 'EconLossSmall', 'Econ Losses - Small')
amies$label

amies.r = 
  amies %>% 
  dplyr::slice(2,4,6)

ggplot(amies.r, aes(x = label, y = AMIE)) + 
  geom_pointrange(aes(ymin = `2.5% CI`,
                      ymax = `97.5% CI`),
                  lwd = 1/2, position = 
                    position_dodge(width = 1/2)) + coord_flip() + 
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + 
  theme_ipsum(base_size = 12) + labs(y = 'AMIE', x = 'Interaction Terms') + 
  theme(panel.grid.minor = element_blank()) +
  ggtitle("Figure A4: Average Marginal Interaction Effects")

ggsave(("fig_a4_amie.pdf"), 
       device = cairo_pdf, width = 6, height = 6)
